Distinguishable and Exchangeable Dyads: Multilevel Modelling

Cross-Sectional and Intensive Longitudinal APIM and DIM

Pascal Küng

2025-10-25

Overview

  • Distinguishable Dyads
    • Cross-Sectional APIM
    • Longitudinal APIM
  • Exchangeable Dyads
    • Cross-Sectional APIM
    • Cross-Sectional DIM
    • Equivalence between APIM and DIM
    • Longitudinal DIM
    • Longitudinal APIM

Loading Libraries and setting CmdStan Backend

Some functions used in this presentation are sourced from the files available in the folder 00_R_Functions.

library(tidyverse)
library(brms)
library(bmlm)
library(easystats)
library(DiagrammeR)
library(DHARMa)
library(MASS)
library(purrr)
source(file.path('00_R_Functions', 'ReportModels.R'))
source(file.path('00_R_Functions', 'PrettyTables.R'))
source(file.path('00_R_Functions', 'PrepareData.R'))

options(
  brms.backend = 'cmdstan',
  brms.file_refit = 'on_change'
)

Distinguishable Dyads

Distinguishable Dyads - Cross-Sectional APIM

Distinguishable Dyads - Data: Simulated Dyads

print_df(head(df))
userID coupleID gender communication satisfaction
1_1 1 1 4.851107 4.841332
1_2 1 2 6.004533 5.577165
2_1 2 1 5.881321 7.248741
2_2 2 2 6.433723 6.960293
3_1 3 1 4.283971 6.928699
3_2 3 2 5.516060 6.077688

Distinguishable Dyads - Cross-Sectional APIM

Distinguishable APIM (e.g., Kenny et al., 2006; Kenny & Cook, 1999; Kenny & Kashy, 2011)

Distinguishable Dyads - Preparing Data

df_apim <- df %>%
  # Reshaping to Actor-Partner format (4-field)
  reshape_dyadic_data(
    person_id = 'userID',
    dyad_id = 'coupleID',
    vars_to_reshape = 'communication'
  )  %>%
  mutate(
    # Optional: grand-mean centering
    communication_actor_gmc = communication_actor - mean(communication_actor),
    communication_partner_gmc = communication_partner - mean(communication_partner),
    
    # Create Dummy-Variables for male and female
    is_female = ifelse(gender == 1, 1, 0),
    is_male = ifelse(gender == 2, 1, 0),
  ) 

print_df(head(df_apim))
userID coupleID gender is_male is_female communication satisfaction communication_actor communication_partner communication_actor_gmc communication_partner_gmc
1_1 1 1 0 1 4.851107 4.841332 4.851107 6.004533 -0.1553604 0.9980657
1_2 1 2 1 0 6.004533 5.577165 6.004533 4.851107 0.9980657 -0.1553604
2_1 2 1 0 1 5.881321 7.248741 5.881321 6.433723 0.8748537 1.4272563
2_2 2 2 1 0 6.433723 6.960293 6.433723 5.881321 1.4272563 0.8748537
3_1 3 1 0 1 4.283971 6.928699 4.283971 5.516060 -0.7224960 0.5095933
3_2 3 2 1 0 5.516060 6.077688 5.516060 4.283971 0.5095933 -0.7224960

Distinguishable Dyads - Fitting the Model in BRMS

Model non-independence either via a shared random intercept or correlated residuals (not both). Posterior predictive checks, residual diagnostics and leave one out cross validation (e.g., via loo_compare()) may help inform which one to use.

For other families like ordinal regression with brms, one might work better than the other.

Distinguishable Dyads - Fitting the Model in BRMS

# Option A: Random-intercept model for non-independence

formula <- bf(
  satisfaction ~ 
    
    # Remove global intercept, introduce male and female intercepts.
    0 + is_male + is_female +
    
    # Actor effect for male and female
    communication_actor_gmc:is_male + 
    communication_actor_gmc:is_female + 
    
    # Partner effect for male and female
    communication_partner_gmc:is_male + 
    communication_partner_gmc:is_female +
    
    # Option 0: Random intercept for male and female (correlated)
    # NOT identifiable in the cross-sectional case:
    # (0 + is_male + is_female | coupleID)
    
    # Option A: Shared dyad effect - interpretable dyad-level random intercept. 
    (1 | coupleID), 
  
  # Two distinct residual variances for males vs females:
  sigma = ~ 0 + is_male + is_female
)

priors <- c(
  # prior(normal(2, 3), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "b", dpar = "sigma"),
  prior(student_t(3, 0, 2.5), class = "sd")
)

model_dist_apim_a <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_a') # Cache the model
)

Distinguishable Dyads - Fitting the Model in BRMS

# Option B: Residual CS for non-independence (no random intercept)

formula <- bf(
  satisfaction ~ 
    
    # Remove global intercept, introduce male and female intercepts.
    0 + is_male + is_female +
    
    # Actor effect for male and female
    communication_actor_gmc:is_male + 
    communication_actor_gmc:is_female + 
    
    # Partner effect for male and female
    communication_partner_gmc:is_male + 
    communication_partner_gmc:is_female +
  
    # Option B: Model non-independence by using compound symmetry of the residuals. 
    # Aligns with the SEM conceptualization / literature
    cosy(gr = coupleID),
  
  # Two distinct residual variances for males vs females:
  sigma = ~ 0 + is_male + is_female
)

priors <- c(
  # prior(normal(2, 3), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "b", dpar = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_dist_apim_b <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_b') # Cache the model
)

Distinguishable Dyads - Check Model Convergence and Fit

Check Rhats and Effective Sample Sizes (ESS_tail and ESS_bulk) directly from the brms summary. Additionally you can (using Model B as an example):

rstan::check_hmc_diagnostics(model_dist_apim_b$fit)

Divergences:
0 of 4000 iterations ended with a divergence.

Tree depth:
0 of 4000 iterations saturated the maximum tree depth of 10.

Energy:
E-BFMI indicated no pathological behavior.
loo::pareto_k_table(loo(model_dist_apim_b))

All Pareto k estimates are good (k < 0.7).

Distinguishable Dyads - Check Model Convergence and Fit

plot(model_dist_apim_b, ask = FALSE)

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, 'dens_overlay_grouped', group = 'is_male')

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, 'ecdf_overlay_grouped', group = 'is_male')

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, type = "loo_pit_overlay")

Distinguishable Dyads - Check Model Convergence and Fit

# Custom function to make DHARMa work with brms (see file 'Functions')
DHARMa.check_brms(model_dist_apim_b)

Distinguishable Dyads - Check Model Convergence and Fit

In case of non-convergence, bad residuals or low predictive accuracy, the model is likely miss-specified.

You can try different families that match your data by using the family() argument. For instance, you can easily conduct an ordinal regression by simply setting family = cumulative() and adjusting your priors. For dichotomous outcomes you can use family = bernoulli(). Brms supports a variety of families and link functions for these families. In general, you should use the default link function (check ?brms::family.brmsfit()).

Distinguishable Dyads - Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
is_male 4.60 [4.47, 4.74] 1.001 5163 3375
is_female 5.61 [5.50, 5.72] 1.000 4706 3152
is_male:communication_actor_gmc 1.84 [1.75, 1.94] 1.002 3380 3262
is_female:communication_actor_gmc 1.59 [1.51, 1.67] 1.002 3812 2935
is_male:communication_partner_gmc 0.17 [0.07, 0.27] 1.000 4128 3151
is_female:communication_partner_gmc 0.31 [0.23, 0.39] 1.000 3663 3475
Residual Structure
cosy 0.50 [0.44, 0.56] 1.001 4076 3169
sigma_is_male 0.47 [0.41, 0.53] 1.000 4362 2672
sigma_is_female 0.26 [0.20, 0.32] 1.001 4595 3241

Note that sigma is reported on the log scale (see ?brms::family.brmsfit()). Using exp(value) retrieves the simulated SD values well.

Distinguishable Dyads - (Intensive) Longitudinal APIM

Distinguishable Dyads - L-APIM - Simulated Dataset

userID coupleID diaryday gender is_female provided_support mvpa provided_support_actor provided_support_partner
31_1 31 1 1 1 3.56 51.49 3.56 2.96
31_1 31 2 1 1 2.72 35.53 2.72 3.9
31_1 31 3 1 1 5 79.79 5 5
... ... ... ... ... ... ... ... ...
31_2 31 1 2 0 2.96 43.78 2.96 3.56
31_2 31 2 2 0 3.9 63.6 3.9 2.72
31_2 31 3 2 0 5 86.78 5 5

Distinguishable Dyads - L-APIM - Centering

With repeated measures we need to center variables in order to not conflate levels of analysis:

  • Between-Person: The person-mean in relation to the grand mean
  • Within-Person: Daily fluctuations of an individual from their person-mean

While it may seem like we have 3-levels (couple/person/day), by including the means of both partners and their correlations in the model, all information about the couple-level (level 3) is included. This will become clear when comparing the exchangeable APIM to the DIM. We thus use a 2-level model with the dyad as the level of analysis.

Distinguishable Dyads - L-APIM - Centering

df_long_apim <- bmlm::isolate(
  df_long_apim, 
  by = 'userID', # NOT coupleID
  value = c('provided_support_actor', 'provided_support_partner'),
  which = 'both' 
) %>%
# renaming to avoid confusion with between- and within COUPLE centred DIM variables
  rename(
    provided_support_actor_cwp = provided_support_actor_cw,
    provided_support_partner_cwp = provided_support_partner_cw,
    provided_support_actor_cbp = provided_support_actor_cb,
    provided_support_partner_cbp = provided_support_partner_cb,
  ) %>%
  # Centering time
  mutate(
    diaryday_c = scale(diaryday, center = TRUE, scale = FALSE)
  )

Distinguishable Dyads - L-APIM - Model nlme

nlme model of this data for reference (simplified random effects structure to achieve convergence).

library(nlme)
# nlme model for reference
model1 <- lme(
  fixed = 
    # Intercepts
    mvpa ~ 0 + is_male + is_female +
    
    # Time Slopes
    is_male:diaryday_c +
    is_female:diaryday_c +
    
    # Between-Person APIM
    is_male:provided_support_actor_cbp + 
    is_male:provided_support_partner_cbp +
    
    is_female:provided_support_actor_cbp + 
    is_female:provided_support_partner_cbp +
  
    # Within-Person APIM
    is_male:provided_support_actor_cwp + 
    is_male:provided_support_partner_cwp +
    
    is_female:provided_support_actor_cwp + 
    is_female:provided_support_partner_cwp,
  
  random = 
    ~ 0 + is_male + is_female 
       
    + is_male:diaryday_c + is_female:diaryday_c
   
  # in nlme we need to remove many random slopes to achieve convergence
  # alternatively, we could try to remove correlations between the slopes. 
    #+ is_male:provided_support_actor_cwp 
    #+ is_male:provided_support_partner_cwp
    #+ is_female:provided_support_actor_cwp 
    #+ is_female:provided_support_partner_cwp 
  | coupleID, 
  
  weights = varIdent(form = ~ 1 | gender), # heterogeneous residual variances
  corr = corCompSymm(form = ~ 1 | coupleID/diaryday), # compound symmetry (partner's values on each day)
  
  data = df_long_apim, 
  control = list(maxIter=1000)
)


summary(model1)
Linear mixed-effects model fit by REML
  Data: df_long_apim 
       AIC     BIC    logLik
  37560.52 37718.9 -18755.26

Random effects:
 Formula: ~0 + is_male + is_female + is_male:diaryday_c + is_female:diaryday_c | coupleID
 Structure: General positive-definite, Log-Cholesky parametrization
                     StdDev     Corr                
is_male               8.5191329 is_mal is_fml is_m:_
is_female            10.5608508 -0.030              
is_male:diaryday_c    0.1218733  0.156  0.071       
is_female:diaryday_c  0.2168773 -0.075 -0.310 -0.033
Residual             21.0281595                     

Correlation Structure: Compound symmetry
 Formula: ~1 | coupleID/diaryday 
 Parameter estimate(s):
      Rho 
0.1462979 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | gender 
 Parameter estimates:
       1        2 
1.000000 1.003816 
Fixed effects:  mvpa ~ 0 + is_male + is_female + is_male:diaryday_c + is_female:diaryday_c +      is_male:provided_support_actor_cbp + is_male:provided_support_partner_cbp +      is_female:provided_support_actor_cbp + is_female:provided_support_partner_cbp +      is_male:provided_support_actor_cwp + is_male:provided_support_partner_cwp +      is_female:provided_support_actor_cwp + is_female:provided_support_partner_cwp 
                                          Value Std.Error   DF  t-value p-value
is_male                                57.64354 1.4629562 4131 39.40210  0.0000
is_female                              54.64020 1.7806967 4131 30.68473  0.0000
is_male:diaryday_c                     -0.04526 0.0351885 4131 -1.28632  0.1984
is_female:diaryday_c                   -0.01491 0.0455931 4131 -0.32706  0.7436
is_male:provided_support_actor_cbp      2.71504 2.0489686 4131  1.32508  0.1852
is_male:provided_support_partner_cbp   -2.48759 1.9517338 4131 -1.27456  0.2025
is_female:provided_support_actor_cbp   -2.90163 2.3204813 4131 -1.25044  0.2112
is_female:provided_support_partner_cbp -1.21943 2.4362413 4131 -0.50054  0.6167
is_male:provided_support_actor_cwp      7.26209 0.4956258 4131 14.65236  0.0000
is_male:provided_support_partner_cwp   -0.63969 0.4999640 4131 -1.27947  0.2008
is_female:provided_support_actor_cwp    6.85210 0.5021507 4131 13.64550  0.0000
is_female:provided_support_partner_cwp -0.61715 0.4964572 4131 -1.24310  0.2139
 Correlation: 
                                       is_mal is_fml is_m:_ is_f:_
is_female                              -0.016                     
is_male:diaryday_c                      0.083  0.038              
is_female:diaryday_c                   -0.055 -0.230  0.063       
is_male:provided_support_actor_cbp     -0.059  0.002  0.000  0.000
is_male:provided_support_partner_cbp    0.055 -0.002  0.000  0.000
is_female:provided_support_actor_cbp   -0.002  0.054  0.000  0.000
is_female:provided_support_partner_cbp  0.002 -0.058  0.000  0.000
is_male:provided_support_actor_cwp      0.000  0.000  0.030  0.003
is_male:provided_support_partner_cwp    0.000  0.000  0.012  0.001
is_female:provided_support_actor_cwp    0.000 -0.001  0.002  0.009
is_female:provided_support_partner_cwp  0.000  0.000  0.004  0.023
                                       is_ml:prvdd_spprt_ctr_cb
is_female                                                      
is_male:diaryday_c                                             
is_female:diaryday_c                                           
is_male:provided_support_actor_cbp                             
is_male:provided_support_partner_cbp    0.179                  
is_female:provided_support_actor_cbp   -0.006                  
is_female:provided_support_partner_cbp -0.034                  
is_male:provided_support_actor_cwp      0.001                  
is_male:provided_support_partner_cwp   -0.006                  
is_female:provided_support_actor_cwp    0.003                  
is_female:provided_support_partner_cwp -0.001                  
                                       is_ml:prvdd_spprt_prtnr_cb
is_female                                                        
is_male:diaryday_c                                               
is_female:diaryday_c                                             
is_male:provided_support_actor_cbp                               
is_male:provided_support_partner_cbp                             
is_female:provided_support_actor_cbp   -0.034                    
is_female:provided_support_partner_cbp -0.006                    
is_male:provided_support_actor_cwp      0.003                    
is_male:provided_support_partner_cwp   -0.001                    
is_female:provided_support_actor_cwp    0.000                    
is_female:provided_support_partner_cwp -0.001                    
                                       is_fml:prvdd_spprt_ctr_cb
is_female                                                       
is_male:diaryday_c                                              
is_female:diaryday_c                                            
is_male:provided_support_actor_cbp                              
is_male:provided_support_partner_cbp                            
is_female:provided_support_actor_cbp                            
is_female:provided_support_partner_cbp  0.179                   
is_male:provided_support_actor_cwp      0.001                   
is_male:provided_support_partner_cwp    0.000                   
is_female:provided_support_actor_cwp    0.001                   
is_female:provided_support_partner_cwp -0.007                   
                                       is_fml:prvdd_spprt_prtnr_cb
is_female                                                         
is_male:diaryday_c                                                
is_female:diaryday_c                                              
is_male:provided_support_actor_cbp                                
is_male:provided_support_partner_cbp                              
is_female:provided_support_actor_cbp                              
is_female:provided_support_partner_cbp                            
is_male:provided_support_actor_cwp      0.000                     
is_male:provided_support_partner_cwp   -0.002                     
is_female:provided_support_actor_cwp    0.014                     
is_female:provided_support_partner_cwp -0.003                     
                                       is_ml:prvdd_spprt_ctr_cw
is_female                                                      
is_male:diaryday_c                                             
is_female:diaryday_c                                           
is_male:provided_support_actor_cbp                             
is_male:provided_support_partner_cbp                           
is_female:provided_support_actor_cbp                           
is_female:provided_support_partner_cbp                         
is_male:provided_support_actor_cwp                             
is_male:provided_support_partner_cwp   -0.084                  
is_female:provided_support_actor_cwp   -0.012                  
is_female:provided_support_partner_cwp  0.145                  
                                       is_ml:prvdd_spprt_prtnr_cw
is_female                                                        
is_male:diaryday_c                                               
is_female:diaryday_c                                             
is_male:provided_support_actor_cbp                               
is_male:provided_support_partner_cbp                             
is_female:provided_support_actor_cbp                             
is_female:provided_support_partner_cbp                           
is_male:provided_support_actor_cwp                               
is_male:provided_support_partner_cwp                             
is_female:provided_support_actor_cwp    0.144                    
is_female:provided_support_partner_cwp -0.012                    
                                       is_fml:prvdd_spprt_ctr_cw
is_female                                                       
is_male:diaryday_c                                              
is_female:diaryday_c                                            
is_male:provided_support_actor_cbp                              
is_male:provided_support_partner_cbp                            
is_female:provided_support_actor_cbp                            
is_female:provided_support_partner_cbp                          
is_male:provided_support_actor_cwp                              
is_male:provided_support_partner_cwp                            
is_female:provided_support_actor_cwp                            
is_female:provided_support_partner_cwp -0.086                   

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.552859205 -0.676127064  0.004364928  0.676509893  3.747382582 

Number of Observations: 4180
Number of Groups: 38 

Distinguishable Dyads - L-APIM - Model brms

Equivalent brms model with simplified random effects structure.

formula <- bf(
  mvpa ~ 0 + is_male + is_female +
    
    # Time Slopes
    is_male:diaryday_c +
    is_female:diaryday_c +
    
    # Between-Person APIM
    is_male:provided_support_actor_cbp + 
    is_male:provided_support_partner_cbp +
    
    is_female:provided_support_actor_cbp + 
    is_female:provided_support_partner_cbp +
  
    # Within-Person APIM
    is_male:provided_support_actor_cwp + 
    is_male:provided_support_partner_cwp +
    
    is_female:provided_support_actor_cwp + 
    is_female:provided_support_partner_cwp +
  
  # Accounting for non-independence between partners' means and trajectories 
  # and effect sensitivities via random effects:
  (0 + is_male + is_female
    + is_male:diaryday_c + is_female:diaryday_c
  | coupleID) +
    
    
  # Accounting for daily non-independence (mimicking nlme's corCompSym)
  # modelled via common 'shocks' with a random intercept:
  (1 | coupleID:diaryday)

  # heteroscedastic residuals
  , sigma ~ 0 + is_male + is_female 
)

priors <- c(
  prior(normal(50, 20), class = "b", coef = "is_male"), # male intercept
  prior(normal(50, 20), class = "b", coef = "is_female"), # female intercept
  prior(normal(0, 5), class = "b") # other beta coefficients
  # Other priors can be specified, but brms choosed defaults well. 
  # Defaults can be seen via default_prior(formula, data)
)

model_dist_apim_long_simple <- brm(
  formula = formula, 
  data = df_long_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_long_simpel') # Cache the model
)

summary(model_dist_apim_long_simple)

Distinguishable Dyads - L-APIM - Results (scrollable)

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
is_male 57.60 [54.56, 60.72] 1.003 799 1522
is_female 54.59 [50.90, 58.43] 1.002 896 1501
is_male:diaryday_c -0.05 [-0.11, 0.02] 1.002 4136 3086
is_female:diaryday_c -0.01 [-0.11, 0.08] 1.003 2332 2377
is_male:provided_support_actor_cbp 2.34 [-1.48, 6.14] 1.004 1131 1857
is_male:provided_support_partner_cbp -2.22 [-5.80, 1.77] 1.005 1025 2022
is_female:provided_support_actor_cbp -2.27 [-6.38, 2.16] 1.005 1131 2003
is_female:provided_support_partner_cbp -0.91 [-5.27, 3.39] 1.004 1173 1369
is_male:provided_support_actor_cwp 7.18 [ 6.23, 8.08] 1.001 5137 2898
is_male:provided_support_partner_cwp -0.68 [-1.65, 0.30] 1.000 5478 3149
is_female:provided_support_actor_cwp 6.80 [ 5.85, 7.77] 1.003 5255 3076
is_female:provided_support_partner_cwp -0.63 [-1.57, 0.33] 1.002 5192 3021
Random Effects
coupleID: sd(is_male) 8.92 [ 6.92, 11.78] 1.004 1135 1880
coupleID: sd(is_female) 10.95 [ 8.70, 14.48] 1.003 1296 1908
coupleID: sd(is_male:diaryday_c) 0.08 [ 0.00, 0.19] 1.006 607 1754
coupleID: sd(is_female:diaryday_c) 0.22 [ 0.13, 0.33] 1.002 1597 1948
coupleID: cor(is_male,is_female) -0.03 [-0.36, 0.30] 1.002 829 1388
coupleID: cor(is_male,is_male:diaryday_c) 0.14 [-0.60, 0.76] 1.000 4203 2417
coupleID: cor(is_female,is_male:diaryday_c) 0.08 [-0.61, 0.73] 0.999 5075 2549
coupleID: cor(is_male,is_female:diaryday_c) -0.06 [-0.45, 0.35] 1.003 2419 2812
coupleID: cor(is_female,is_female:diaryday_c) -0.27 [-0.63, 0.14] 1.001 2480 2262
coupleID: cor(is_male:diaryday_c,is_female:diaryday_c) -0.03 [-0.76, 0.72] 1.023 217 592
coupleID:diaryday: sd(Intercept) 8.02 [ 6.74, 9.20] 1.007 701 1424
Residual Structure
sigma_is_male 2.97 [ 2.94, 3.01] 1.001 1604 2549
sigma_is_female 2.97 [ 2.93, 3.01] 1.004 1618 2618

Distinguishable Dyads - L-APIM - AR1 Model

We can try to add an AR1 residual structure.

    ... +
    ar(time = diaryday, gr = coupleID:userID, p = 1),
    
  sigma ~ 0 + is_male + is_female
)

Distinguishable Dyads - L-APIM - AR1 Model

AR1 can changes the fixed effects quite a bit. You may want to check which model performs better! Often the same day shock term and the time slopes are enough.

loo1 <- loo(model_dist_apim_long_simple)
loo2 <- loo(model_dist_apim_long_ar)

loo::pareto_k_table(loo1)

All Pareto k estimates are good (k < 0.7).
loo::pareto_k_table(loo2)

All Pareto k estimates are good (k < 0.7).

Distinguishable Dyads - L-APIM - AR1 Model

a <- loo_compare(loo1, loo2)
print(a)
                            elpd_diff se_diff
model_dist_apim_long_simple    0.0       0.0 
model_dist_apim_long_ar     -107.9      25.3 
report::report(a)
The difference in predictive accuracy, as indexed by Expected Log Predictive
Density (ELPD-LOO), suggests that 'model_dist_apim_long_simple' is the best
model (ELPD = -18689.18), followed by 'model_dist_apim_long_ar' (diff-ELPD =
-107.87 +- 25.27, p < .001)

In this instance, the AR1 model performs much worse! We keep the simpler model.

Maximal Model Specification

We add all random slopes and their correlations. This is not converging in nlme.

formula <- bf(
  mvpa ~ 
    
    # Intercepts
    0 + is_male + is_female +
    
    # Time Slopes
    is_male:diaryday_c +
    is_female:diaryday_c +
    
    # Between-Person APIM
    is_male:provided_support_actor_cbp + 
    is_male:provided_support_partner_cbp +
    
    is_female:provided_support_actor_cbp + 
    is_female:provided_support_partner_cbp +
  
    # Within-Person APIM
    is_male:provided_support_actor_cwp + 
    is_male:provided_support_partner_cwp +
    
    is_female:provided_support_actor_cwp + 
    is_female:provided_support_partner_cwp +
    
    # Accounting for non-independence between partners' means and trajectories 
    # and effect sensitivities via random effects:
    (0 + is_male + is_female +  
       
       is_male:diaryday_c + is_female:diaryday_c +
       
       is_male:provided_support_actor_cwp + 
       is_male:provided_support_partner_cwp +
       is_female:provided_support_actor_cwp + 
       is_female:provided_support_partner_cwp 
    | coupleID ) + 
    
    
    # Accounting for daily non-independence
    (1 | coupleID:diaryday), 
  
  # No more AR1. You could also test MA or ARMA. 
    
  sigma ~ 0 + is_male + is_female   # heteroscedastic residuals
)

priors <- c(
  prior(normal(50, 20), class = "b", coef = "is_male"),
  prior(normal(50, 20), class = "b", coef = "is_female"),
  prior(normal(0, 5), class = "b")
)


model_dist_apim_long_complex <- brm(
  formula = formula, 
  data = df_long_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_long_complex') # Cache the model
)

Maximal Model Specification: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
is_male 57.69 [54.55, 60.79] 1.006 670 1461
is_female 54.64 [50.74, 58.38] 1.003 855 1253
is_male:diaryday_c -0.04 [-0.11, 0.02] 1.001 3824 2732
is_female:diaryday_c -0.01 [-0.10, 0.08] 1.001 2527 2714
is_male:provided_support_actor_cbp 2.35 [-1.64, 6.14] 1.002 1150 1631
is_male:provided_support_partner_cbp -2.19 [-6.08, 1.76] 1.002 1079 1766
is_female:provided_support_actor_cbp -2.26 [-6.75, 2.09] 1.007 1124 1916
is_female:provided_support_partner_cbp -1.00 [-5.28, 3.33] 1.001 1155 1944
is_male:provided_support_actor_cwp 7.17 [ 6.17, 8.18] 1.000 5204 3280
is_male:provided_support_partner_cwp -0.69 [-1.84, 0.40] 1.000 4397 3071
is_female:provided_support_actor_cwp 6.75 [ 5.76, 7.76] 1.000 5223 3141
is_female:provided_support_partner_cwp -0.61 [-1.65, 0.46] 1.001 4853 3519
Random Effects
coupleID: sd(is_male) 9.00 [ 6.94, 11.96] 1.002 1236 2065
coupleID: sd(is_female) 11.08 [ 8.68, 14.67] 1.006 1082 1715
coupleID: sd(is_male:diaryday_c) 0.10 [ 0.00, 0.20] 1.006 723 1365
coupleID: sd(is_female:diaryday_c) 0.22 [ 0.14, 0.33] 1.002 1444 2487
coupleID: sd(is_male:provided_support_actor_cwp) 0.55 [ 0.02, 1.88] 1.002 1946 1865
coupleID: sd(is_male:provided_support_partner_cwp) 1.14 [ 0.06, 2.80] 1.003 1268 1488
coupleID: sd(is_female:provided_support_actor_cwp) 0.56 [ 0.02, 1.90] 1.001 1863 1850
coupleID: sd(is_female:provided_support_partner_cwp) 0.79 [ 0.04, 2.45] 1.001 1388 1989
coupleID: cor(is_male,is_female) -0.03 [-0.35, 0.29] 1.006 684 1184
coupleID: cor(is_male,is_male:diaryday_c) 0.08 [-0.44, 0.60] 1.000 4681 2884
coupleID: cor(is_female,is_male:diaryday_c) 0.05 [-0.48, 0.55] 1.002 4532 2350
coupleID: cor(is_male,is_female:diaryday_c) -0.05 [-0.43, 0.32] 1.001 2947 2950
coupleID: cor(is_female,is_female:diaryday_c) -0.23 [-0.56, 0.16] 1.001 2767 2782
coupleID: cor(is_male:diaryday_c,is_female:diaryday_c) -0.03 [-0.61, 0.55] 1.005 612 1123
coupleID: cor(is_male,is_male:provided_support_actor_cwp) 0.08 [-0.56, 0.67] 1.002 6756 2939
coupleID: cor(is_female,is_male:provided_support_actor_cwp) 0.03 [-0.62, 0.63] 1.000 5597 2577
coupleID: cor(is_male:diaryday_c,is_male:provided_support_actor_cwp) -0.02 [-0.65, 0.61] 1.001 4911 3277
coupleID: cor(is_female:diaryday_c,is_male:provided_support_actor_cwp) -0.04 [-0.64, 0.61] 1.001 5031 2968
coupleID: cor(is_male,is_male:provided_support_partner_cwp) 0.06 [-0.53, 0.60] 1.000 5161 2838
coupleID: cor(is_female,is_male:provided_support_partner_cwp) 0.14 [-0.47, 0.64] 1.002 4251 2529
coupleID: cor(is_male:diaryday_c,is_male:provided_support_partner_cwp) 0.23 [-0.46, 0.76] 1.002 2351 3107
coupleID: cor(is_female:diaryday_c,is_male:provided_support_partner_cwp) -0.12 [-0.66, 0.53] 1.001 4085 3266
coupleID: cor(is_male:provided_support_actor_cwp,is_male:provided_support_partner_cwp) 0.01 [-0.63, 0.62] 1.000 3047 3703
coupleID: cor(is_male,is_female:provided_support_actor_cwp) 0.05 [-0.61, 0.67] 1.001 6151 2751
coupleID: cor(is_female,is_female:provided_support_actor_cwp) 0.15 [-0.52, 0.69] 1.001 5501 2835
coupleID: cor(is_male:diaryday_c,is_female:provided_support_actor_cwp) 0.04 [-0.60, 0.65] 1.000 4526 2979
coupleID: cor(is_female:diaryday_c,is_female:provided_support_actor_cwp) -0.05 [-0.65, 0.59] 1.001 5504 3576
coupleID: cor(is_male:provided_support_actor_cwp,is_female:provided_support_actor_cwp) 0.03 [-0.59, 0.65] 1.001 3388 3383
coupleID: cor(is_male:provided_support_partner_cwp,is_female:provided_support_actor_cwp) 0.07 [-0.59, 0.67] 1.000 3143 2769
coupleID: cor(is_male,is_female:provided_support_partner_cwp) 0.08 [-0.52, 0.62] 1.002 5913 3145
coupleID: cor(is_female,is_female:provided_support_partner_cwp) -0.03 [-0.59, 0.57] 1.001 5643 3114
coupleID: cor(is_male:diaryday_c,is_female:provided_support_partner_cwp) -0.07 [-0.68, 0.57] 1.001 3847 3483
coupleID: cor(is_female:diaryday_c,is_female:provided_support_partner_cwp) 0.07 [-0.57, 0.65] 1.000 4974 2750
coupleID: cor(is_male:provided_support_actor_cwp,is_female:provided_support_partner_cwp) 0.06 [-0.58, 0.70] 1.000 2694 3363
coupleID: cor(is_male:provided_support_partner_cwp,is_female:provided_support_partner_cwp) 0.01 [-0.59, 0.64] 1.000 3436 3302
coupleID: cor(is_female:provided_support_actor_cwp,is_female:provided_support_partner_cwp) -0.01 [-0.64, 0.64] 1.001 2498 3420
coupleID:diaryday: sd(Intercept) 8.00 [ 6.60, 9.18] 1.010 583 1433
Residual Structure
sigma_is_male 2.97 [ 2.93, 3.01] 1.003 1404 2008
sigma_is_female 2.97 [ 2.93, 3.01] 1.004 1327 2137

Maximal Model Specification: Diagnostics

In complex models like this, we want to make sure we are not overfitting.

  • Conduct regular model convergence checks (inspecting chains, Rhats, ESS, Multicollinearity, Residual diagnostics etc.). If these are bad (but good for a simpler model), it could be a sign of overfitting.

  • Leave-one-out cross-validation is a powerful tool for investigating overfitting. It tests out of sample prediction and thus punishes overfitting.

Maximal Model Specification: Diagnostics

loo1 <- loo(model_dist_apim_long_simple)
loo2 <- loo(model_dist_apim_long_complex)

loo::pareto_k_table(loo1)

All Pareto k estimates are good (k < 0.7).
loo::pareto_k_table(loo2)
Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     4179  100.0%  299     
   (0.7, 1]   (bad)         1    0.0%  <NA>    
   (1, Inf)   (very bad)    0    0.0%  <NA>    

High Pareto-Ks can be a sign of overfitting or other problems. Here, both models are fine.

The loo values (elpd) themselves cannot be interpreted individually, but we can compare models (next page)

Maximal Model Specification: Diagnostics

a <- loo_compare(loo1, loo2)
print(a)
                             elpd_diff se_diff
model_dist_apim_long_simple   0.0       0.0   
model_dist_apim_long_complex -3.9       1.6   

The difference between the models seems quite small. No obvious overfitting, but also no improvement over the simpler model (even slightly worse, but hard to say for such small differences). Some rule of thumbs say that if the difference is higher than 2x the SE or 4x the SE we can say there is a difference.

In our case, we could follow our theory and the philosophy of the maximal random effects structure (Barr et al., 2013). Alternatively, we could follow the philosophy of parsimony.

(continued on next page)

Maximal Model Specification: Diagnostics

One approach is to use the report package’s approach to obtain a p-value for the comparison (see their documentation for which assumptions they make):

report::report(a)
The difference in predictive accuracy, as indexed by Expected Log Predictive
Density (ELPD-LOO), suggests that 'model_dist_apim_long_simple' is the best
model (ELPD = -18689.18), followed by 'model_dist_apim_long_complex' (diff-ELPD
= -3.87 +- 1.59, p = 0.015)

In case of clear overfitting, we could try to improve the model by subsequently excluding the smallest random effects (or correlations between them) and compare various models. For more guidance on random effects check out del Rosario & West (2025).

Exchangeable Dyads

Exchangeable Dyads - Cross-Sectional APIM

Exchangeable APIM (e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Cross-Sectional APIM

Exchangeable APIM (e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Main Assumptions

Partners are exchangeable, i.e., not systematically different.

  • Equal actor effects
  • Equal partner effects
  • Equal means
  • Equal residual variances

But they should still be allowed to vary within each couple, while being correlated.

(e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Cross-Sectional APIM: Data

Same cross-sectional data in the same 4-field actor-partner format as before.

print_df(head(df_apim))
userID coupleID gender is_male is_female communication satisfaction communication_actor communication_partner communication_actor_gmc communication_partner_gmc
1_1 1 1 0 1 4.851107 4.841332 4.851107 6.004533 -0.1553604 0.9980657
1_2 1 2 1 0 6.004533 5.577165 6.004533 4.851107 0.9980657 -0.1553604
2_1 2 1 0 1 5.881321 7.248741 5.881321 6.433723 0.8748537 1.4272563
2_2 2 2 1 0 6.433723 6.960293 6.433723 5.881321 1.4272563 0.8748537
3_1 3 1 0 1 4.283971 6.928699 4.283971 5.516060 -0.7224960 0.5095933
3_2 3 2 1 0 5.516060 6.077688 5.516060 4.283971 0.5095933 -0.7224960

Exchangeable Dyads - Cross-Sectional APIM: Model

formula <- bf(
  satisfaction ~ 1 + 
    communication_actor_gmc + communication_partner_gmc + 
    
    # Option 1: A single couple level random intercept.
    # In the cross sectional case this is the maximal random effects structure and
    # we cannot distinguish each partners variance from noise. 
    # (1 | coupleID)
    
    # Option 2: using a compound symmetry residual structure. 
    cosy(gr = coupleID)
  
  # Note: no need to model separate sigmas for each partner.
  # Homogeneous residual variance is estimated. 
  # Implied: sigma = ~ 1
)

priors <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_ind_apim <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_apim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional APIM: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 5.11 [5.00, 5.22] 1.000 4801 3093
communication_actor_gmc 1.72 [1.66, 1.78] 1.000 5048 3183
communication_partner_gmc 0.24 [0.18, 0.30] 1.001 5120 2749
Residual Structure
cosy 0.34 [0.27, 0.41] 1.002 4332 2705
sigma 1.55 [1.49, 1.62] 1.002 4215 2828

Exchangeable Dyads - Cross-Sectional APIM: Test for Distinguishability

Leave-one-out (loo) cross-validation for model comparison (even if not nested).

a <- loo_compare(
  loo(model_ind_apim), 
  loo(model_dist_apim_b)
)
print(a)
Model elpd_diff se_diff
model_dist_apim 0.0 0.0
model_ind_apim -172.9 16.5
report::report(a)

The difference in predictive accuracy, as indexed by Expected Log Predictive Density (ELPD-LOO), suggests that ‘model_dist_apim_b’ is the best model (ELPD = -1806.42), followed by ‘model_ind_apim’ (diff-ELPD = -172.85 +- 16.53, p < .001).
See: documentation of “report”

Exchangeable Dyads - Cross-Sectional DIM

Exchangeable Dyads - Cross-Sectional DIM: Data

Starting from scratch (no 4-field data needed)

userID coupleID communication satisfaction
1_1 1 4.851107 4.841332
1_2 1 6.004533 5.577165
2_1 2 5.881321 7.248741
2_2 2 6.433723 6.960293
3_1 3 4.283971 6.928699
3_2 3 5.516060 6.077688

Exchangeable Dyads - Cross-Sectional DIM: Centering

Decompose communication variance into:

  • Between-couple (cbc): Couple-mean communication skills in relation to other couples
  • Within-couple (cwc): Individuals’ communication skills in relation to their couple-mean

Same assumption about exchangeability as in the APIM

Exchangeable Dyads - Cross-Sectional DIM: Centering

df_dim <- df_apim2 %>%
  group_by(coupleID) %>%
  mutate(
    communication_cm = mean(communication, na.rm = TRUE),
    communication_cwc= communication - communication_cm
  ) %>%
  ungroup() %>%
  mutate(
    communication_cbc= communication_cm - mean(communication_cm, na.rm = TRUE)
  ) %>% 
  dplyr::select(c('userID', 'coupleID', 'satisfaction', 'communication_cwc', 'communication_cbc'))

print_df(head(df_dim))
userID coupleID satisfaction communication_cwc communication_cbc
1_1 1 4.841332 -0.5767130 0.4213527
1_2 1 5.577165 0.5767130 0.4213527
2_1 2 7.248741 -0.2762013 1.1510550
2_2 2 6.960293 0.2762013 1.1510550
3_1 3 6.928699 -0.6160447 -0.1064514
3_2 3 6.077688 0.6160447 -0.1064514

Exchangeable Dyads - Cross-Sectional DIM: Model

formula <- bf(
  satisfaction ~ 1 + 
    communication_cbc+ communication_cwc + 
    cosy(gr = coupleID) 
)

priors <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_ind_dim <- brm(
  formula = formula, 
  data = df_dim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_dim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional DIM: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 5.11 [5.01, 5.21] 1.003 5216 2981
communication_cbc 1.96 [1.88, 2.04] 1.000 4664 3077
communication_cwc 1.48 [1.39, 1.57] 1.003 4715 3276
Residual Structure
cosy 0.34 [0.26, 0.41] 1.000 4312 3208
sigma 1.55 [1.48, 1.62] 1.000 4479 3342

Equivalence APIM and DIM

APIM (left) vs. DIM (right)

\[b_{actor\_gmc} + b_{partner\_gmc} = b_{cbc}\] \[b_{actor\_gmc} - b_{partner\_gmc} = b_{cwc}\]

(Bolger et al., 2025)

Equivalence APIM and DIM

Results

Top sliders: grand-mean-centered Communication for Actor and Partner (APIM).
Bottom sliders: reparameterized communication — Centered Between-Couple (xcbc) and Centered Within-Couple (xcwc) (DIM).

0.00 0.00 0.00 0.00

Equivalence APIM and DIM

  • If the couple mean goes up by 1 and the within-couple deviation from the couple mean stays fixed, this means that both partners’ predictors must go up by 1 unit. Thus, the effects are linear combinations:

\[b_{cbc} = b_{actor\_gmc} + b_{partner\_gmc}\]

  • Similarly, if a person’s deviation from their couple mean is one unit higher and the couple mean stays fixed, this means their partners’ deviation must be one unit lower. Thus, the effects are a linear combination:

\[b_{cwc} = b_{actor\_gmc} - b_{partner\_gmc}\]

(Bolger et al., 2025)

Equivalence APIM and DIM

  • Conversely, to obtain the actor and partner effects given the DIM estimates:

\[b_{actor\_gmc} = \frac{b_{cbc} + b_{cwc}}{2}\]

\[b_{partner\_gmc} = \frac{b_{cbc} - b_{cwc}}{2}\]

Equivalence APIM and DIM

Example: using APIM coefficients to obtain DIM coefficients and computing Credible Intervals of DIM coefficients.

a <- hypothesis(
  model_ind_apim, 
  "communication_actor_gmc + communication_partner_gmc = 0"
)

a$hypothesis[,c(2,3,4,5)]
  Estimate  Est.Error CI.Lower CI.Upper
1 1.961219 0.04314024 1.876668 2.044644

APIM (left) vs. DIM (right)

Equivalence APIM and DIM: Takeaway

  • APIM and DIM are reparametrizations of the same model
  • APIM: intuitive actor/partner framing
  • DIM: clean separation of between vs within components
  • Estimating one model allows for directly obtaining estimates of the other
  • Same random-effects structure at dyad level and same assumptions

Side Note:

In the distinguishable case, the Dyadic Score Model (DSM) (e.g., Stadler et al., 2023) is equivalent to the APIM! (Bolger et al., 2025; Iida et al., 2018)

Exchangeable Dyads - Longitudinal APIM/DIM: Simulated Data

coupleID userID diaryday diaryday_c mvpa provided_support_actor_cwp provided_support_partner_cwp provided_support_actor_cbp provided_support_partner_cbp
31 31_1 1 -27 51.49 -0.1 -0.44 0.75 0.49
31 31_1 2 -26 35.53 -0.95 0.5 0.75 0.49
31 31_1 3 -25 79.79 1.33 1.6 0.75 0.49
... ... ... ... ... ... ... ... ...
31 31_2 1 -27 43.78 -0.44 -0.1 0.49 0.75
31 31_2 2 -26 63.6 0.5 -0.95 0.49 0.75
31 31_2 3 -25 86.78 1.6 1.33 0.49 0.75

Exchangeable Dyads - Longitudinal DIM/APIM

  • We need an APIM or DIM on both the between person level and the within-person level.
  • Due to equivalence, we could have a DIM on one level and an APIM on the other.

If we want a within-person DIM:

df_long_dim <- df_long_apim2 %>%
  mutate(
    
    # Between Person Level DIM (centering between couples and within couple between person)
    provided_support_cbc = (provided_support_actor_cbp + provided_support_partner_cbp) / 2,
    provided_support_cwcbp = (provided_support_actor_cbp - provided_support_partner_cbp) / 2,
    
    # Within Person Level DIM (couple mean and deviation on each day)
    provided_support_cwp_mean = (provided_support_actor_cwp + provided_support_partner_cwp) / 2,
    provided_support_cwp_halfdiff = (provided_support_actor_cwp - provided_support_partner_cwp) / 2
  ) 

Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the sum and difference approach

We need to randomly assign -1 and +1 for each Partner within each couple. This will be needed for contrast coding the intercept.

For example:

df_long_dim <- df_long_dim %>%
  group_by(coupleID) %>%
  mutate(
    base = ifelse(userID == min(userID), 1, -1),
    flip = 1 - 2 * rbinom(1, 1, 0.5),   # yields +1 or -1 once per couple
    Idiff = base * flip
  ) %>%
  ungroup() %>%
  relocate(Idiff, .after = userID) %>%
  dplyr::select(-base, -flip)


print_df(print_couple_preview(df_long_dim, couple_id = "31"))  

Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the Sum and Difference Approach

coupleID userID Idiff diaryday diaryday_c mvpa provided_support_actor_cwp provided_support_partner_cwp provided_support_actor_cbp provided_support_partner_cbp provided_support_cbc provided_support_cwcbp provided_support_cwp_mean provided_support_cwp_halfdiff
31 31_1 -1 1 -27 51.49 -0.1 -0.44 0.75 0.49 0.62 0.13 -0.27 0.17
31 31_1 -1 2 -26 35.53 -0.95 0.5 0.75 0.49 0.62 0.13 -0.22 -0.73
31 31_1 -1 3 -25 79.79 1.33 1.6 0.75 0.49 0.62 0.13 1.47 -0.13
31 31_1 -1 4 -24 61.52 -0.18 1.6 0.75 0.49 0.62 0.13 0.71 -0.89
... ... ... ... ... ... ... ... ... ... ... ... ... ...
31 31_2 1 1 -27 43.78 -0.44 -0.1 0.49 0.75 0.62 -0.13 -0.27 -0.17
31 31_2 1 2 -26 63.6 0.5 -0.95 0.49 0.75 0.62 -0.13 -0.22 0.73
31 31_2 1 3 -25 86.78 1.6 1.33 0.49 0.75 0.62 -0.13 1.47 0.13
31 31_2 1 4 -24 91.18 1.6 -0.18 0.49 0.75 0.62 -0.13 0.71 0.89

Exchangeable Dyads - Longitudinal APIM/DIM: Model

formula <- bf(
  mvpa ~ 1 + 
    
    diaryday_c +
  
    # Within-person APIM
    provided_support_actor_cwp + 
    provided_support_partner_cwp +
    # Equivalent to within-person DIM:
    # provided_support_cwp_mean + 
    # provided_support_cwp_halfdiff +
  
    # Between-person APIM
    provided_support_actor_cbp +
    provided_support_partner_cbp +
    # Equivalent to between-person DIM
    # provided_support_cwcbp+
    # provided_support_cbc +
    
    # Dyad-Level intercept and slopes for time-varying predictors
    (1 + diaryday_c + provided_support_actor_cwp + provided_support_partner_cwp | coupleID) +
    # Both partners' deviations from these dyad-level means and slopes.
    # Note: separate random effects block to make them uncorrelated form dyad-level REs. 
    (0 + Idiff + 
         I(Idiff * diaryday_c) + 
         I(Idiff * provided_support_actor_cwp) + 
         I(Idiff * provided_support_partner_cwp) | coupleID) +
    
    # Accounting for same-day shocks/coupling
    (1 | coupleID:diaryday)
    
    # Autocorrelated residuals
    # Due to bad fit in the distinguishable case, we omit AR1 from the start.
    # ar(time = diaryday, gr = coupleID:userID, p = 1)
  
  # Again, no need to model heterogeneous residual variances (sigma)
  # Implied: sigma = ~ 1
)

priors <- c(
  prior(normal(50, 20), class = "Intercept"),
  prior(normal(0, 10), class = "b")
  
  # We could set priors for other things, but brms sets good default priors
  #prior(student_t(3, 0, 1.5), class = "sd"),
  #prior(lkj(2), class = "cor"), # correlation prior for random effect matrix
  #prior(student_t(3, 0, 1.5), class = "sigma")
)

model_apim_ind_long <- brm(
  formula = formula, 
  data = df_long_dim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'model_apim_ind_long_apim')
)

Exchangeable Dyads - Longitudinal APIM/DIM: Explanation Idiff

For appropriate random effects, we can use the Sum and Difference approach (del Rosario & West, 2025; Kenny & Ackerman, 2023):

  • Randomly give one partner a constant of -1 and the other a constant of 1
  • The couple level intercept represents the mean (or sum) of both partners’ intercepts. (1 | coupleID)
  • A column of 1s and -1s represent deviations (difference) of each partner from the couple intercept with each partner contributing equally in one direction (+1 and -1)

Exchangeability is retained: if partners are flipped, the result is the same.

(del Rosario & West (2025) provide practical guidance on how to reduce the random effects structure in case of non-convergence.)

Exchangeable Dyads - Longitudinal APIM/DIM: Results

Comparing APIM (left) and DIM (right) output

APIM (left) vs. DIM (right) APIM

Equivalence holds on both levels. Numeric deviations arise due to priors, sampling noise, and uncertainty in some of the effects in this example.

Exchangeable Dyads - Longitudinal APIM/DIM: Results

APIM Results in Detail:

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 56.37 [53.90, 58.65] 1.004 721 1478
diaryday_c -0.03 [-0.09, 0.03] 1.001 3224 2931
provided_support_actor_cwp 7.03 [ 6.30, 7.78] 1.000 4444 3193
provided_support_partner_cwp -0.63 [-1.41, 0.16] 0.999 4486 2972
provided_support_actor_cbp 0.11 [-3.12, 3.28] 1.001 746 1596
provided_support_partner_cbp -2.16 [-5.36, 1.19] 1.003 833 1130
Random Effects
coupleID: sd(Intercept) 7.14 [ 5.54, 9.40] 1.003 1040 1678
coupleID: sd(diaryday_c) 0.11 [ 0.02, 0.19] 1.001 818 714
coupleID: sd(provided_support_actor_cwp) 0.59 [ 0.03, 1.68] 1.003 1272 2016
coupleID: sd(provided_support_partner_cwp) 0.77 [ 0.04, 1.98] 1.003 1040 1508
coupleID: sd(Idiff) 7.26 [ 5.69, 9.49] 1.003 835 1509
coupleID: sd(IIdiffMUdiaryday_c) 0.13 [ 0.07, 0.19] 1.000 1536 1410
coupleID: sd(IIdiffMUprovided_support_actor_cwp) 0.36 [ 0.02, 1.14] 1.001 2480 2465
coupleID: sd(IIdiffMUprovided_support_partner_cwp) 0.84 [ 0.05, 1.87] 1.000 1372 1696
coupleID: cor(Intercept,diaryday_c) -0.10 [-0.57, 0.40] 1.000 2353 1828
coupleID: cor(Intercept,provided_support_actor_cwp) 0.27 [-0.62, 0.85] 1.001 4069 2694
coupleID: cor(diaryday_c,provided_support_actor_cwp) -0.07 [-0.81, 0.77] 1.000 3677 2612
coupleID: cor(Intercept,provided_support_partner_cwp) 0.24 [-0.58, 0.83] 1.002 3836 2599
coupleID: cor(diaryday_c,provided_support_partner_cwp) 0.08 [-0.72, 0.80] 1.002 2664 2519
coupleID: cor(provided_support_actor_cwp,provided_support_partner_cwp) 0.18 [-0.76, 0.86] 1.002 2091 2733
coupleID: cor(Idiff,IIdiffMUdiaryday_c) -0.13 [-0.53, 0.30] 1.002 3067 2395
coupleID: cor(Idiff,IIdiffMUprovided_support_actor_cwp) 0.19 [-0.71, 0.86] 1.000 5906 2490
coupleID: cor(IIdiffMUdiaryday_c,IIdiffMUprovided_support_actor_cwp) -0.12 [-0.84, 0.75] 1.000 4577 3129
coupleID: cor(Idiff,IIdiffMUprovided_support_partner_cwp) -0.20 [-0.80, 0.58] 1.001 4121 2846
coupleID: cor(IIdiffMUdiaryday_c,IIdiffMUprovided_support_partner_cwp) 0.48 [-0.45, 0.92] 1.002 2829 2575
coupleID: cor(IIdiffMUprovided_support_actor_cwp,IIdiffMUprovided_support_partner_cwp) -0.15 [-0.85, 0.77] 1.003 2287 2860
coupleID:diaryday: sd(Intercept) 8.00 [ 6.62, 9.15] 1.003 683 1185
Residual Structure
sigma 19.49 [18.88, 20.14] 1.001 1015 1972

Extract Full APIM Random Effects Variance-Covariance Matrix

We can rotate the random effects structure back to obtain the full APIM Random effects variance-covariance matrix (Kenny & Ackerman, 2023). Just as we would in a SEM model with equality constraints!

For example, the within-person variance of the intercept is:

\[var(Intercept) + var(Idiff)\]

and the cross-person covariance of both partners’ intercepts is:

\[var(Intercept) - var(Idiff)\]

Extract Full APIM Random Effects Variance-Covariance Matrix

# Custom function 
rotate_apim_covariance <- function(
    fit,
    gr = "coupleID",
    Idiff = "Idiff"
) {
  random_effects <- fit$ranef
  random_effects <- random_effects[random_effects$group == gr,]

  SUM <- random_effects$coef[!grepl(Idiff, random_effects$coef)]
  DIFF <- random_effects$coef[grepl(Idiff, random_effects$coef)]
  
  varcor <- VarCorr(fit)[[gr]]
  
  within_person_covariance_matrix <- varcor$cov[SUM,'Estimate',SUM] + varcor$cov[DIFF,'Estimate',DIFF]
  cross_person_covariance_matrix <- varcor$cov[SUM,'Estimate',SUM] - varcor$cov[DIFF,'Estimate',DIFF]
  
  p <- nrow(within_person_covariance_matrix)
  
  Full <- rbind(
    cbind(within_person_covariance_matrix,  cross_person_covariance_matrix),
    cbind(cross_person_covariance_matrix, within_person_covariance_matrix)
  )
  
  # nice labels
  base <- SUM
  rn <- c(paste0("PartnerA_", base), paste0("PartnerB_", base))
  colnames(Full) <- rn
  rownames(Full) <- rn
  
  return(list(within_person_covariance_matrix=within_person_covariance_matrix, cross_person_covariance_matrix=cross_person_covariance_matrix, full_covariance_matrix=Full))
}

a <- rotate_apim_covariance(model_apim_ind_long)

Extract Full APIM Random Effects Variance-Covariance Matrix: Within-Person Matrix

print_df(a$within_person_covariance_matrix)
Intercept diaryday_c provided_support_actor_cwp provided_support_partner_cwp
Intercept 108.0688313 -0.1972175 2.0422566 0.0893135
diaryday_c -0.1972175 0.0310245 -0.0126549 0.0612627
provided_support_actor_cwp 2.0422566 -0.0126549 0.9063907 0.0523327
provided_support_partner_cwp 0.0893135 0.0612627 0.0523327 1.9091184

Extract Full APIM Random Effects Variance-Covariance Matrix: Cross-Person Matrix

print_df(a$cross_person_covariance_matrix)
Intercept diaryday_c provided_support_actor_cwp provided_support_partner_cwp
Intercept -2.1620176 0.0421946 0.6281852 2.7491930
diaryday_c 0.0421946 -0.0028062 0.0014472 -0.0434522
provided_support_actor_cwp 0.6281852 0.0014472 0.3676112 0.2182846
provided_support_partner_cwp 2.7491930 -0.0434522 0.2182846 0.0067028

Extract Full APIM Random Effects Variance-Covariance Matrix: Full Matrix

print_df(a$full_covariance_matrix)
PartnerA_Intercept PartnerA_diaryday_c PartnerA_provided_support_actor_cwp PartnerA_provided_support_partner_cwp PartnerB_Intercept PartnerB_diaryday_c PartnerB_provided_support_actor_cwp PartnerB_provided_support_partner_cwp
PartnerA_Intercept 108.0688313 -0.1972175 2.0422566 0.0893135 -2.1620176 0.0421946 0.6281852 2.7491930
PartnerA_diaryday_c -0.1972175 0.0310245 -0.0126549 0.0612627 0.0421946 -0.0028062 0.0014472 -0.0434522
PartnerA_provided_support_actor_cwp 2.0422566 -0.0126549 0.9063907 0.0523327 0.6281852 0.0014472 0.3676112 0.2182846
PartnerA_provided_support_partner_cwp 0.0893135 0.0612627 0.0523327 1.9091184 2.7491930 -0.0434522 0.2182846 0.0067028
PartnerB_Intercept -2.1620176 0.0421946 0.6281852 2.7491930 108.0688313 -0.1972175 2.0422566 0.0893135
PartnerB_diaryday_c 0.0421946 -0.0028062 0.0014472 -0.0434522 -0.1972175 0.0310245 -0.0126549 0.0612627
PartnerB_provided_support_actor_cwp 0.6281852 0.0014472 0.3676112 0.2182846 2.0422566 -0.0126549 0.9063907 0.0523327
PartnerB_provided_support_partner_cwp 2.7491930 -0.0434522 0.2182846 0.0067028 0.0893135 0.0612627 0.0523327 1.9091184

Convert to correlation matrix of random effects

# Custom function
sd_cor <- function(Sigma) {
  sds  <- sqrt(diag(Sigma))
  cors <- cov2cor(Sigma)  
  list(sd = sds, cor = cors)
}

full <- sd_cor(a$full_covariance_matrix)

Random effects SDs

print_df(as.data.frame(round(full$sd, 3)))
round(full$sd, 3)
PartnerA_Intercept 10.396
PartnerA_diaryday_c 0.176
PartnerA_provided_support_actor_cwp 0.952
PartnerA_provided_support_partner_cwp 1.382
PartnerB_Intercept 10.396
PartnerB_diaryday_c 0.176
PartnerB_provided_support_actor_cwp 0.952
PartnerB_provided_support_partner_cwp 1.382

Random effect correlation matrix

print_df(as.data.frame(round(full$cor, 3)))
PartnerA_Intercept PartnerA_diaryday_c PartnerA_provided_support_actor_cwp PartnerA_provided_support_partner_cwp PartnerB_Intercept PartnerB_diaryday_c PartnerB_provided_support_actor_cwp PartnerB_provided_support_partner_cwp
PartnerA_Intercept 1.000 -0.108 0.206 0.006 -0.020 0.023 0.063 0.191
PartnerA_diaryday_c -0.108 1.000 -0.075 0.252 0.023 -0.090 0.009 -0.179
PartnerA_provided_support_actor_cwp 0.206 -0.075 1.000 0.040 0.063 0.009 0.406 0.166
PartnerA_provided_support_partner_cwp 0.006 0.252 0.040 1.000 0.191 -0.179 0.166 0.004
PartnerB_Intercept -0.020 0.023 0.063 0.191 1.000 -0.108 0.206 0.006
PartnerB_diaryday_c 0.023 -0.090 0.009 -0.179 -0.108 1.000 -0.075 0.252
PartnerB_provided_support_actor_cwp 0.063 0.009 0.406 0.166 0.206 -0.075 1.000 0.040
PartnerB_provided_support_partner_cwp 0.191 -0.179 0.166 0.004 0.006 0.252 0.040 1.000

References

Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278. https://doi.org/10.1016/j.jml.2012.11.001
Bolger, N., Laurenceau, J.-P., & DiGiovanni, A. (2025). Unified analysis model for indistinguishable and distinguishable dyads. Innovations in Interpersonal Relationships and Health Research: Advancing the Integration of Interdisciplinary Approaches to Dyadic Behavior Change. https://doi.org/10.17605/OSF.IO/WYDCJ
Bürkner, P.-C. (2024). Brms: Bayesian regression models using stan. https://github.com/paul-buerkner/brms
del Rosario, K. S., & West, T. V. (2025). A Practical Guide to Specifying Random Effects in Longitudinal Dyadic Multilevel Modeling. Advances in Methods and Practices in Psychological Science, 8(3), 25152459251351286. https://doi.org/10.1177/25152459251351286
Gabry, J., Češnovar, R., Johnson, A., & Bronder, S. (2024). Cmdstanr: R interface to CmdStan. https://mc-stan.org/cmdstanr/
Guo, J., Gabry, J., Goodrich, B., Johnson, A., Weber, S., & Badr, H. S. (2025). Rstan: R interface to stan. https://mc-stan.org/rstan/
Hartig, F. (2024). DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. http://florianhartig.github.io/DHARMa/
Iida, M., Seidman, G., & Shrout, P. E. (2018). Models of interdependent individuals versus dyadic processes in relationship research. Journal of Social and Personal Relationships, 35(1), 59–88. https://doi.org/10.1177/0265407517725407
Kenny, D. A., & Ackerman, R. A. (2023). Estimation of random effects with over-time dyadic data using multilevel modeling: The sum and difference method. OSF Preprints. https://doi.org/10.31219/osf.io/fju72
Kenny, D. A., & Cook, W. (1999). Partner effects in relationship research: Conceptual issues, analytic difficulties, and illustrations. Personal Relationships, 6(4), 433–448. https://doi.org/10.1111/j.1475-6811.1999.tb00202.x
Kenny, D. A., & Kashy, D. A. (2011). Dyadic Data Analysis Using Multilevel Modeling. In Handbook of Advanced Multilevel Analysis. Routledge.
Kenny, D. A., Kashy, D. A., & Cook, W. L. (2006). Dyadic data analysis. Guilford Press.
Lüdecke, D., Makowski, D., Ben-Shachar, M. S., Patil, I., Wiernik, B. M., Bacher, E., & Thériault, R. (2025). Easystats: Framework for easy statistical modeling, visualization, and reporting. https://easystats.github.io/easystats/
Stadler, G., Scholz, U., Bolger, N., Shrout, P. E., Knoll, N., & Lüscher, J. (2023). How is companionship related to romantic partners’ affect, relationship satisfaction, and health behavior? Using a longitudinal dyadic score model to understand daily and couple-level effects of a dyadic predictor. Applied Psychology: Health and Well-Being, 15(4), 1530–1554. https://doi.org/10.1111/aphw.12450
Vehtari, A., Gabry, J., Magnusson, M., Yao, Y., Bürkner, P.-C., Paananen, T., & Gelman, A. (2024). Loo: Efficient leave-one-out cross-validation and WAIC for bayesian models. https://mc-stan.org/loo/
Vuorre, M. (2023). Bmlm: Bayesian multilevel mediation. https://github.com/mvuorre/bmlm/
Wickham, H. (2023). Tidyverse: Easily install and load the tidyverse. https://tidyverse.tidyverse.org